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(~| Abstract. We consider the motion of uncharged dust grains of arbitrary shape including 
I ' the effects of electromagnetic radiation and thermal emission. The resulting relativisti- 

O ; 

^ I cally covariant equation of motion is expressed in terms of standard optical parameters. 
Explicit expressions for secular changes of osculating orbital elements are derived in de- 



tail for the special case of the Poynting-Robertson effect. Two subcases are considered: 
(i) central acceleration due to gravity and the radial component of radiation pressure 
independent of the particle velocity, (ii) central acceleration given by gravity and the 
radiation force as the disturbing force. The latter case yields results which may be com- 
pared with secular orbital evolution in terms of orbital elements for an arbitrarily shaped 
dust particle. The effects of solar wind are also presented. 
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1. Introduction 

Poynting (1904) formulated a problem of finding equation of motion for a perfectly ab- 
sorbing spherical particle under the action of electromagnetic radiation. Poynting did not 
succeed in finding correct solution. The second case, closely connected with the rclativis- 
tic equation of motion for a free particle under the action of electromagnetic radiation, 
was presented by Einstein (1905), who calculated the change of energy and light pres- 
sure at arbitrary angle of incidence on a plane mirror. Robertson (1937) subsequently 
derived a correct equation of motion for a perfectly absorbing spherical particle. This 
result has been applied to astronomical problems for several decades and is known as 
the Foynting-Robertson (P-R) effect. Robertson's case was relativistically generalized by 
Klacka (1992a), who showed that the generalized P-R effect holds only for the special 
case where the total momentum of the outgoing radiation per unit time is colinear with 
the incident radiation (in the proper/rest frame of reference of the particle), which may 
include radiation normally incident on a plane mirror. Since real particles interact with 
electromagnetic radiation in a complicated manner and particles of various optical prop- 
erties exist (e. g., Mishchenko et al. 2002), it is essential to have an equation of motion 
sufficiently general to cover a wide range of optical parameters, not just the limited cases 
previously investigated. As an attempt of formulating such a general equation of motion, 
we can mention Lyttlcton (1976), Klacka (1994a), Klacka and Kocifaj (1994), Kocifaj et 
al. (2000). The last three presentations are of hypothetical character only, since none of 
them proves correctness of the equation of motion. However, Klacka (2000a) succeeded 
in deriving relativistically covariant equation of motion. This equation of motion was 
derived in other possible ways by Klacka (2000b), Klacka and Kocifaj (2001), within 
the accuracy corresponding to Klacka (1994a). Finally, knowing the results obtained by 
Klacka and Kocifaj (1994) and Kocifaj et al. (2000) and having in disposal papers by 
Klacka (2000a, 2000b, 2000c), Kimura et al. (2002) presented a repetition of the equation 
of motion corresponding to Klacka and Kocifaj (1994) and Kocifaj et al. (2000). 

The equations of motion for a moving particle have been constructed under the as- 
sumption that thermal emission from a particle is isotropic and does not exert a radiation 
pressure force on the particle in the particle frame of reference. Recently, Mishchenko 
(2001) has formulated the radiation pressure on arbitrarily shaped particles arising from 
an anisotropy of thermal emission. In this paper, we derive the equation of motion for an 
arbitrarily shaped particle moving in a radiation field taking into account the radiation 
pressure caused by an anisotropy of thermal emission as well as scattering and absorp- 
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tion of light. Relativistically covariant equation of motion is presented. Derivation of the 
equation of motion is physically fully reasoned. 

We begin by reviewing in Sec. 2 and 3 the basic physical processes in proper and 
stationary frames. The equation of motion for simultaneous action of gravity and elec- 
tromagnetic radiation is presented in Sec. 4. Sec. 5 then applies these results to the 
calculation of osculating orbital elements, including the special case of the P-R effect, 
which may be treated analytically. The calculation is carried out in Sec. 6 to first order 
in v/c and applied to the ejection of a dust particle from a parent body such as a comet 
or asteroid. Two cases are considered: (i) the disturbing acceleration is given in terms of 
velocity (Robertson 1937, Wyatt and Whipple 1950), and (ii) the electromagnetic radi- 
ation itself is a disturbing function. However, these authors (including Lyttleton 1976) 
did not obtain the correct expression for the secular change of longitude of pericenter 
(perihelion). In Sec. 7 we use the equation of motion to second order in vjc for the P-R 
effect. In particular we obtain a correct expression for the secular change of longitude of 
pericenter (perihelion). Sec. 8 then finds the secular change in the advancement of perihe- 
lion to first order in w/c, with gravity as a central acceleration. Next, we briefly discuss, 
in Sec. 9, the secular evolution of orbital elements for the P-R effect and nearly-circular 
orbits. Sec. 10 then treats the effect of the solar wind on the secular changes in the orbital 
elements to the second order in v ju^ where u is the speed of solar wind particles. Finally, 
Sec. 11 summarizes our results. 

Other theoretical papers on the basic properties of the P-R effect were written during 
the last decades: Burns et al. (1979), Mediavilla and Buitrago (1989), Mignard (1992), 
Srikanth (1999), Williams (2002). Some others will be mentioned within the context of 
the discussed problems. Since some confusion exists in presenting derivations and results 
in the most cited papers (Robertson 1937; Wyatt and Whipple 1950; Burns et al. 1979; 
Mignard 1992), we present detailed derivation of the secular changes of orbital elements 
for heliocentric orbits and the P-R effect. 

2. Proper reference freime of the particle — stationary particle 

The term "stationary particle" will denote a particle which does not move in a given 
inertial frame of reference. Primed quantities will denote quantities measured in the 
proper reference frame of the particle rest frame of the particle. 

The flux density of photons scattered into an elementary solid angle dfl' = 
sinO' d9' d(f)' is proportional to p'{9',(j>') dft' , where p'{6',(f)') is the "phase function". 
The phase function depends on orientation of the particle with respect to the direction 
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of the incident radiation and on the particle characteristics; angles 9' , (/>' correspond to 
the direction (and orientation) of travel of the scattered radiation, 6' is the polar angle 

which vanishes for propagation along the unit vector e'l of the incident radiation. The 
phase function fulfills the normalisation condition 

/ p'{6',(i>')dn' = 1. (1) 

Jin 

The momentum of the incident beam of photons which is lost in the process of inter- 
action with the particle is proportional to the cross-section C^^^ (extinction). The part 
proportional to C^^^, (absorption) is emitted in the form of thermal radiation and the 
part proportional to C'^^i — C^^^, = C^^a i'' scattered. The differential scattering cross 
section dC'^^^/d^' = C^^^ p'{d' , cf)') depends on the polarization state of the incident light 
as well as on the incidence and scattering directions (e. g., Mishchenko et al. 2002). 

The momentum (per unit time) of the scattered photons into an elementary solid 
angle dO.' is 

dp' sea = -S' CL, v'{0', <!>') K' dfl' , (2) 
c 

where the unit vector in the direction of scattering is 

K' = cose' e[ + sine' cos(t)' e'^ + sine' sine/)' eg . (3) 

>S" is the flux density of radiation energy (energy flow through unit area perpendicular to 
the ray per unit time). The system of unit vectors used on the RHS of the last equation 
forms an orthogonal basis. The total momentum (per unit time) of the scattered photons 
is 

P'.ca = - 5' C^ca / p'{e',(l>')K' dn' . (4) 

The momentum (per unit time) obtained by the particle due to the interaction with 
radiation - radiation force acting on the particle - is 



(5) 



^ = i 5' {CL* e; - C7L„ £ p'{e',<l>') K' rffi'l + F',{T') , 

where the emission component of the radiation force acting on the particle of absolute 
temperature T' is (Mishchenko et al. 2002, pp. 63-66) 

F',iT') = -- r div' f r' K{r',T',J) dr' . (6) 

The unit vector r = r' jr' is given by position vector r' of the observation point with 
origin inside the particle (the emitted radiation in the far-field zone of the particle prop- 
agates in the radial direction, i. e., along the unit vector f'), w' is (angular) frequency of 
radiation, 
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K'^{r',T',u,')=lUT',iv') 



K[,{r',uj')- [ Z[,{r',r",u;')dr" 

Ji-K 



(7) 



/ 47r 

where K'^-^ is the (1,1) element of the particle extinction matrix, Z'^^ is the (1,1) element 
of the phase matrix and the Planck blackbody energy distribution is given by the well- 
known relation 

'i<->'>™{-(ii;)-f • 

Thermal emission has to be included in the total interaction of the particle with electro- 
magnetic radiation: if the particle's absolute temperature is above zero, it can emit as 
well as scatter and absorb electromagnetic radiation. 

For the sake of brevity, we will use dimensionless efficiency factors instead of cross 
sections C^: C'^ = Q'^ A', where A' is geometrical cross section of a sphere of volume 
equal to the volume of the particle. Equation (5) can be rewritten to the form 



^ = 1 y A' {[Q'^^, - < cose' > Q',J e[ + [- < smO' cos<j>' > Q',J e', + 
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[ - < Sine' sin,/.' > g;,J ea + ^ i^V , (9) 

where < x' >= J^^ x' p'{9' , </)') dW and = F'^{T') ■ e^. As for the energy, we assume 
that it is conserved: the energy (per unit time) of the incoming radiation E[, equals to 
the energy (per unit time) of the outgoing radiation (after interaction with the particle) 
E'g. We will use the fact that time t' = r, where r is proper time. 

Summarizing important equations, we can write them in a short form 

where = Q'^^^ - < cos6i' > Q'2 = - < sinO' coscj)' > Q',^^, Q'^ = - < 

sin 9' smcj)' > Q's^a- We have added an assumption of equilibrium state when the parti- 
cle's mass does not change. 

2.1. Summary of the important equations 

Using the text concerning energy below Eq. (9) and the last Eq. (10), we may describe the 
total process of interaction in the form of the following equations (energies and momenta 
per unit time): 

K = E[ = A' 5' , 

3 

p'o^ {I - g'l) - (Q^ + Q'^ e'^) E'Jc - ^ F',^ e'^ , 

p\ = {E'Jc) e'l , (11) 
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The index "i" represents the incoming radiation, beam of photons, the index "o" rep- 
resents the outgoing radiation. The relation for p'^ represents a generaUzation of the 

following cases: 

i) Q\ = 1, Q2 = Q3 = Fg'i = F'^^ = = - Robertson (1937), Robertson and Noonan 
(1968), Srikanth (1999); 

ii) Q; arbitrary, Q'^ = Q'^ = F^^ = F',^ = F;,^ = - Klacka (1992a); 

iii) Q[, Q'2, Q3 arbitrary, F^^ = F^^ = K3 = " Klacka (2000a). 

The changes of energy and momentum of the particle due to the interaction with 
electromagnetic radiation are 

— - F, - E^-O, 

%-v\-v'.. (12) 



3. Stationary frame of reference 

By the term "stationary frame of reference" (laboratory frame) we mean a frame of 
reference in which particle moves with a velocity vector v = v(€). The physical quantities 
measured in the stationary frame of reference will be denoted by unprimed symbols. 

Our aim is to derive equation of motion for the particle in the stationary frame of 
reference. Wc will use the fact that we know this equation in the proper frame of reference 
- see Eqs. (11) and (12). 

If we have a four- vector = {A^, A), where A° is its time component and A is its 
spatial component, generalized special Lorentz transformation yields 

A°' = j{A° - V A/c) , 

A' = A + [{'y - l)v A/v"^ - 7 A°/c] v , (13) 
with inverse 

AO = 7 {A°' + V ■ A'/c) , 

A = A' + [{j - l)v A'/v^ + 7 A°'/c] V , (14) 

where 7 = 1/ -y/l — v^/c^ . 

As for four-vectors we immediately introduce four-momentum: 



p'' = (pO,p)^(F/c,p) 



(15) 



3.1. Incoming radiation 

Applying Eqs. (14) and (15) to quantity {E[/c,p'i) (four-momentum per unit time - 
proper time is a scalar quantity) and taking into account also Eqs. (11), we can write 

Ei = E[-i {\ + V e'Jc) , 

Pi = ^ {e'l + [(7 - 1) • e'Jv^ + 7/c] v) . (16) 
Using the fact that = {h v, h ei) for photons, we have 

I/' = V Wi , 

e'l = i- {ei + [(7 - 1) V ei/i;2 - 7/c] v] , (17) 
where 



wi = 7 (1 — V ■ ei/c) . (18) 

Inserting the second of Eqs. (17) into Eqs. (18), one obtains 
Ei = E'. , 

= (l/wi) {E'Jc) ei . (19) 

We have four-vector = {Ei/c,pi) = (l,ei)£',/c = {l/w\,ei/w\) wi Ei/c = 
b1 wi Ei/c. 

For monochromatic radiation the flux density of radiation energy becomes 

S' =n' hv' c; S = nhuc, (20) 

where n and n' are concentrations of photons (photon number densities) in the corre- 
sponding reference frames. We also have continuity equation 

^^.j'' = 0, f = {cn,cne^), (21) 
with current density j'^. Application of Eq. (13) then yields 

n' = win. (22) 
Using Eqs. (17), (20) and (22) we finally obtain 

S' = wlS . (23) 
Eqs. (11), (19) and (23) then together give Ei=wi S A', pi = wi S A' e\/c. 



8 



3.2. Outgoing radiation 

The situation is analogous to that of the preceding subsection. It is only a little more 
algebraically complicated, since radiation may also spread out in directions given by unit 
vectors 62, 63. How can we find transformations for the unit vectors e'2 and 63? The 
crucial point is what physics do these unit vectors describe. The vectors e'2 and 63 can 
be used to describe directions of propagation of radiation scattered by the particle. Thus, 
aberration of light also exists for each of these unit vectors. The relations between 62 and 
62, 63 and 63, are analogous to that presented by the second of Eq. (17): 
1 



e;. = — {ej + [(7 - 1) v-Sj/v^ - 7/c] v] , j = 1,2,3, (24) 



where 

Wj=j{l - VBj/c) , J = 1,2,3. (25) 

It is worth mentioning that vectors {e^ ;i = 1,2,3} form an orthonormal set of vectors, 
and, unit vectors {ej;j = 1,2,3} are not orthogonal unit vectors. 

Applying Eqs. (14) and (15) to the quantity {E'^/c,p'g) (four-momentum per unit time 
- proper time is a scalar quantity), we can write 



TP' 

Po = Po + 



(7 - 1) • pljv^ + ^ 5^ 



V . (26) 



Using also p\ = E[ e'^/c and Eqs. (11), (24), (26), 

Eo = g'l wiEi^ + {1 - Q[) Ei + 

wi Ei {Q'2 + Q3) 7 - wi Ei {Q'2/w2 + Qa/ws) - 

3 

^ej (c/Wj - 7 C) , 

j=l 

/-, n,i \ Ei „, Wi Ei ^ , tt'i El 

Po = (1 - Qi) — ei + Qi — — 1 V - Qj — — (c Bj/wj - 7 t;) - 

1 ^ 

^e, (ce,M-7t^) • (27) 
3.3. Equation of motion 

In analogy with Eqs. (12), we have for the changes of energy and momentum of the 
particle due to the interaction with electromagnetic radiation 

dE 

—, — = Ei — Eo , 
d T 

^=P.-P.. (28) 
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Putting Eqs. (27) into Eqs. (28), using also Pi = (Ei/c) ei, one easily obtains 



3 



dp / „, Wi Ei ^ Tp/\ f ^ 



Eq. (29) may be rewritten in terms of four-vectors: 



(30) 



where p** is four- vector of the particle of mass m 

pi^ = mu>^ , (31) 
four- vector of the world- velocity of the particle is 

= (7 c, 7 w) . (32) 
We have also other four- vectors 

6^^ = (lM,e,-M) , i = 1,2,3. (33) 

It can be easily verified that: 

i) the quantity w is a scalar quantity sec first of Eqs. (19); 

ii) Eq. (30) reduces to Eq. (10) for the case of proper inertial frame of reference of the 
particle; 

iii) Eq. (30) yields d mid r = 0. 
We introduce 

6° = \lw, = 7 (1 + ^ • / c) , 

hj = Bj I Wj = e; + [{^ - l)ve'^ / + J / c]v , j = 1,2,3 (34) 

for the purpose of practical calculations. Physics of these relations corresponds to aber- 
ration of light. 

We have derived an equation of motion for real dust particle under the action of 
electromagnetic radiation (including thermal emission). It is supposed that the equation 
of motion is represented by Eqs. (11) and (12) in the proper frame of reference of the 
particle. The final covariant form is represented by Eq. (30), or using Ei = wi S A' (see 
Eqs. (11), (19) and (23)), 

dp^' fw^ S A' 1 , \ , ^ „, 

= ^ ( + c ~ ■ ^^^^ 
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Another form of equation of motion is presented in Appendix A. There is also ex- 
plained why considerations presented by Kimura et al. (2002) are incorrect. 

To first order in v/c, Eqs. (34)-(35) yield 

/ 3 

= ^ ^'i ii^ - ei/c + V ■ e,/c) e, - v/c] + 

ej = (1 - v-e'j/c)e'j + v/c , j = 1,2,3. (36) 

It is worth mentioning to stress that the values of Q'— coefficients depend on particle's 
orientation with respect to the incident radiation - their values are time dependent. A 
little different derivation of the equation of motion within the accuracy to the first order in 
i;/c is presented in Klacka and Kocifaj (2001). Inspite of knowing the correct theoretical 
results, Kimura et al. (2002) have succeeded in publishing a paper where different set of 
unit vectors is used (in reality, it is repetition of the older results of Klacka and Kocifaj 
1994, Kocifaj et al. 2000). The relation between our set of unit vectors {ej;j = 1,2,3} 
and the set of unit vectors {kj;j = 1,2,3} used by Kimura et al. (2002) is following: 
ei = ki, 62 = k2 + {v- k2/c) {k\ — k2) + v/c, 63 = fes — (t; • fes/c) k^ + v/c and definition 
t; • fes = is used - this transformation yields Eqs. (1) - (3) in Kimura et al. (2002): 
e'i = {1 + v ■ k\/c)k\_ — v/c, e'2 = k2 + {v ■ k2/c)k\, 63 = fes. As for physics of the paper 
Kimura et al. (2002), we refer to Appendix C. 

It can be verified that Eq. (35) (or Eqs. (36) within the accuracy to the first order 
in w/c) yields as special cases the situations discussed in Einstein (1905) and Robertson 
(1937) - Robertson's case is obtained simply by substituting Q'l = 1, Q'2 = Q3 ~ 0, 
= = P'e3 = 0) Einstein's results require a little more calculations (see Appendix 

B). 

3.4.. Heuristic derivation 

Since we completely understand the physics of Eq. (30), we are able to present a short 

simple derivation of Eq. (36). We have: 

i) dv'/dt = {[S'A'/{mc)] Q[ + F'^^/m] e[ (see Eq. (10)); dm/dt = is supposed; 

ii) S' = S{\ ~ 2v ■ ei/c) (sec Eq. (23)), due to the change of concentration of photons 
n' ^ n{l - v ■ ei/c) Eq. (22) and Dopplcr effect ly' ^ - v ■ ei/c) - Eq. (17); 

ui) e'j — {1 + v ■ ej/c)ej — v/c, j = 1, 2, 3 (aberration of fight - see Eq. (24)). 

Taking into account these physical phenomena, we finally obtain Eq. (36). However, only 
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relativistic covariant formulation proves that {Qj,F^j;j = 1,2,3} behave as relativisti- 
cally invariant quantities. 

A reader may compare this heuristic derivation of the general case with heiiristic 
derivation of the special case {Q[ = 1, Q2 = Q3 = 0, F'^ = 0) presented by Burns et al. 
(1979, pp. 5-6). 

Another instructive text enabling better understanding of physics of the P-R effect 
can be found in (http://xxx.lanl.gov,astro-ph/0108210), which deals with explanation 
presented in Harwit (1988; pp. 176-177). 

A little different derivation of the relation between S' and S, based on: a) stress- 
energy tensor (energy- momentum tensor), b) transformations of electric and magnetic 
fields, may be found in Klacka (1992a: sections 2.3 and 2.5). 

3.5. Continuous distribution of density flux of energy 

For a continuous frequency distribution of density flux of energy, we can write 

- i: «; + k) e', . (37) 

Taking into account that concentration of photons fulfills n' = wi n (Eq. (22)) and 
that Doppler effect yields v' = wi v (Eq. (17)), we have dn' /dv' = dnjdv. Lorentz 
transformation finally yields 




- E q;- + \ - ■ (38) 



As a consequence, dm/dr = (this corresponds to the condition dE'/dr = 0). 
If the particle is not irradiated, then one has to use 




instead of Eq. (38). As a consequence, Eq. (39) yields {dm/dT)e — ^ej I Mass 

of the particle decreases due to the thermal emission, alone. 



4. Gravitation and electromagnetic radiation 

The generally covariant equation of motion can be immediately written on the basis of 
Eq. (38) (see e. g., Landau and Lifshitz 1975): 
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i=i 



wl S A' 



Q'. + Fi, {b^ - u^/c) 



(40) 



where the operator D / dr is the "total" covariant derivative in the general relativistic 
sense and includes gravitational effects. 

4-1- Gravity and radiation - equation of motion to the first order in v/c 

We can write, on the basis of Eqs. (40), (25) and (33) 
G Mq 



d V 

~dt 



G M:, 



ei + 



V ■ ei V ■ e,. 



+ 



1 

m 







v 













(41) 



where 
= 

/3i = 

/32 = 

0s = 
w 



I - «-e;/c)e;. + v/c, J = 1,2,3, 

p2 /-oo 

^ / Bo(A) {CLt(A/«;) - C:,„(A/«;) ^UA/w^)} dA 
3 m c 7o 

50(A){ -C;,,(A/H <?^(A/u;)}dA, 



G Mq m. c Jo 

-0 



G M0 m c Jq 
1 — V ■ ei / c , 



BQ{\){-C',,,{\/w)g'^{X/w)}d\, 



(42) 



if we use Sun as a source of gravitation and radiation. Rq denotes the radius of the 
Sun and ^©(A) is the solar radiance at a wavelength of A; G, Mq, and r are the grav- 
itational constant, the mass of the Sun, and the distance of the particle from the cen- 
ter of the Sun, respectively. The asymmetry parameter vector g' is defined by g' = 
(1/Ggg^) / n' {dGg^^/dCl')dfl' , where n' is a unit vector in the direction of scattering; g' 
= g[ e'l + g'2 e'2 + ^3 63, ej = {l + v- Cj/c) ej — v/c, • e'j = Stj. We may mention 
that /3i = /32 = /? Q2/Q1, Pa = P Qs/O'i correspond to quantities used in Klacka and 
Kocifaj (2001). 



5. Motion and orbital elements 

Equation of motion is given by Eq. (40), or, by Eq. (41) to the first order in v/c. We 
have to take into account that the non-dimensional parameter, for Solar System 

LqA' 



G Mq m c 



^, ^ 0.02868 ^, A' [m^] 
AnG Mqtuc^^ 12 TT m [kg] 



(43) 
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may change during the motion. Here Lq is the rate of energy outflow from the Sun, the 
solar luminosity. /3 may change and at each point of the orbit all three parameters (3i, 
f32 and /?3 have to be numerically calculated (except for a very special cases, when the /3 
parameters are constant during the motion). 

If we are interested in orbital evolution in terms of osculating orbital elements, we 
come to the crucial point: Which type of osculating orbital elements is correct? Mainly, if 
the "radial radiation pressure" (the dominant term containing parameter /3 = /3i in Eq. 
(75)) has to be considered together with the central gravitational force, or not (compare 
Klacka 1992b). 

As for a definition of osculating orbital elements, we refer to any textbook on celestial 
mechanics. Brouwer and Clemence (1961) write (p. 273): "As the motion progresses under 
the infiuence of the various attracting bodies, the coordinates and velocity components 
at any instant may be used to obtain a set of six orbital elements. These are precisely the 
elements of the ellipse that the body would follow if from that particular instant on, the 
accelerations caused by all "perturbing" bodies ceased to exist." As for our purposes, it 
is sufficient to make a small change: "attracting bodies" are replaced by "forces" . 

In reality, our perturbing physical force corresponds to the total electromagnetic ra- 
diation force. Thus, physics incites us to use gravitational force alone, as the central ac- 
celeration. However, it may occur that the dominant term containing parameter (3i = [3 
in Eq. (41) is comparable to central gravity term. As a consequence, orbital elements will 
change very rapidly during the motion. This would suggest to divide physical force into 
two parts and to use "radial radiation pressure" as a part of central acceleration, i. e. 
together with central gravitational acceleration. The final problem concerns the changes 
of (3 (almost random changes) which prevent us to use the term "Keplerian orbit" for an 
unperturbed problem. 

Thus, we have to use gravitational acceleration of the central body (Sun) as an ac- 
celeration defining Keplerian orbits - complete electromagnetic radiative acceleration is 
a disturbing acceleration - in Eq. (41). As for long-term orbital evolution, we may use 
mean values of orbital elements - they may be defined as time mean values of orbital 
elements when true anomaly changes in 2 tt radians. 

Let us consider that dust particle is ejected from a parent body. Orbital elements 
of the parent body, at the moment of ejection, are: op, ep, ip, fJp, ojp and Op, i. e., 
semi-major axis, eccentricity, inclination, longitude of the ascending node, longitude of 
pericenter - longitude of perihelion for the case of Solar System - and position angle of 
the body measured from the ascending node in the direction of the motion of the body. 
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The velocity vector of the parent body, in a given coordinate system with the origin in 
the dominant central body (Sun in the case of the Solar System) is: 

vp = vr epR + vt epT , 



I GMq . 
y ap (1 - ej,) 

'"T = \ r^ 2T (1 + cosfp) , 

\j ap {1 ~ ej,) 

ap (l - ej,) 

rp = — J- epR 

1 + ep cos jp 



epR = (cos Op cos Op — sinfip sin6p coszp , 

sinfip cosGp + cos Op sinOp cosip , sin Op sinip) , 
epT = (— cosOp sinOp — sinOp cosGp cosip , 

— sinOp sinOp + cosOp cosOp cosip , cosGp sinip) , 
fp = Qp - wp . (44) 

Radial and transversal velocity components are vr and Vt, unit vector epx is normal 
to the radial vector and oriented in orientation of motion. True anomaly fp equals for 
pericenter /perihelion and n for apocenter/aphelion. 

The initial conditions for the particle ejected with velocity A from the parent body 
are: 

Vin = rp , 

Vin = VP + A . (45) 

Eqs. (44)-(45) immediately yield that initial orbital elements of the particle are equal 
to those of the parent body for A = 0: 

Oi„ = Op , LOin = U)p , Qin = ©P • (46) 

(Initial orbital elements for general case represented by Eq. (45) can be calculated from 
the relations presented in Eq. (47) below.) 

Complete equation of motion for dust particle orbiting the central body of mass Mq 
under central's body electromagnetic radiation and gravity interaction is given by Eq. 
(41). We can find osculating orbital elements for dust particle according to the following 
equations: 
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H = r X V 




e = 1 + 



2pE 
GMq 



P 



q = a {1 — e) 



a = 



1 



i = arccos(-ff2/|iy|) , 
sinCl sini = Hx/\H\ , — cosO sini = Hy/\H\ 
eR = r/\r\ = {x,y,z)/r , 
eN = H/\H\ , 

sin 6 sin i = z/r , 

cos6 sini = {bt)^ = {y - x Hy) / {r\H\) , 



where q is perihelion distance. If i = 0, then 6 has to be obtained from cos(0 + 6) = x/r, 
sin(0 + 6) = y/r, and, we may take O in an arbitrary way (e. g., the requirement that O 
is a continuous function of time may define f2 for i = 0). As for the value of inclination, 
we repeat the well-known fact: prograde orbits exist for i £ (0,7r/2), retrograde orbits 
exist for i G {tt/2, tt). The case e > 1 describes the case when the particle is ejected from 
the Solar System. 

5.1. Poynting- Robertson effect 

Putting Q'^ = Q'^ = F'^^ = F'^^ = Ks = ^ in Eq. (40), we obtain Poynting-Robertson 
effect (Robertson 1937, Klacka 1992a): 



Eq. (48) enables us some analytical calculations for changes of osculating orbital 
elements: it is supposed that Q'^ is a constant. We will present them, in the following 
sections. Section 6 will consider the first order in u/c of Eq. (48), section 7 will consider 
also the second order in v/c of Eq. (48). 




cos(e - w) = V ■ btI (e^/G Mq/p) - 1, 



(47) 




(48) 



6. P-R effect — equation of motion to the first order in v/c 

We can write, on the basis of Eqs. (48), (41) - (43) 
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77 = -^«« + ^/3{(l-^) (49) 

where /j. = G Mq, en = e\ and /3 = /3iis a non-dimensional parameter ("the ratio of 
radiation pressure force to the gravitational force"; see also Eq. (43)). Eq. (43) reduces to 
/3 = 5.7 X lQ~^Q'-^^/{Q[g/cm^] s[cm]), for homogeneous spherical particle: q is mass density 
and s is radius of the sphere. 



6.1. Secular changes of orbital elements - radiation pressure as a part of central acceler- 
ation 

We have to use — /x (1 — /3) ej? / as a central acceleration determining osculating 
orbital elements if we want to take a time average (T is time interval between passages 
through two following pericenters) in an analytical way 



Vm (1 - /3) 1 r 



"0 



3/2 9^ / S'(//3) / dfl3 

ai'^ 27r Jo y/fj, (1 - /?) 



a2 yrr 



e2 27r JO 



27r 



5(//3) H d/^ , (50) 



assuming non-pseudo-circular orbits and the fact that orbital elements exhibit only small 
changes during the time interval T; a^j is semi-major axis, ep is eccentricity, fp is true 
anomaly, pp = a^(l — e^); the second and the third Kepler's laws were used: dfp/dt = 
sj IJ,{1 — (3)pj3 ~ conservation of angular momentum, a^/T"^ = /Lt(l — /3)/(47r^). 
Rewriting Eq. (49) into the form 

dv /u (1 - /3) 



dt 



we can immediately write for components of perturbation acceleration to Keplerian mo- 
tion: 

F,H = -2/34^, F,, = -/j4^, F,. = 0, (52) 
c c 

where Fp r, Fp t and Fp n arc radial, transversal and normal components of perturbation 
acceleration, and two-body problem yields 



m(i-/?) . , 

V Pl3 



vpT = \ r^^ (1 + e^cos/ff) . (53) 

V P0 
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The important fact that perturbation acceleration is proportional to v/c (<C 1) ensures 
the above metioned small changes of orbital elements during the time interval T. 

Perturbation equations of celestial mechanics yield for osculating orbital elements {ap 
- semi-major axis; ep eccentricity; inclination (of the orbital plane to the reference 
frame); fi/j longitude of the ascending node; ujj3 - longitude of pcriccntcr; Qfi is the 
position angle of the particle on the orbit, when measured from the ascending node in 
the direction of the particle's motion, 6/3 = W/j + //j): 

^ = 1^ \/ ^(/- /3) ^^^^^^ smf, + F,T{l + e, cos/,)} , 

^ = \ a\ sm f0 + F0T cos//? + f- \ , 

dt y IJ.{1 - /3) [ [ l + e,cos//3jJ 

= -^===== Fp N cos 9/3 

dQp r ^ sin 

dt ~ ^nil - (3)pii sini/3 

dhJf) _ 1 j p'p / _ ^ + '^Z' '^^^ 

dt 





- /3)P/3 
r 




- 








- 

r 








- /3)P/3 



[77 f 77 2 + e/3COs//3 \ 

<^ J> ii cos //3 - F/3 T— f- sm //3 ^ - 

I l + e/3C0s//3 "^J 



„ sin 9/3 

^> TV . . cos 1/3 

smz/3 



c;B/3 - /3)P/3 r- sin 9/3 . 

— = -^^^ 5 - , Fr n —. — - cos^/3 , (54) 

dt r2 - P)Pf3 sini/3 ^ ^ 

where r = pp/{l + ep cos//?). 

Inserting Eqs. (52) (53) into Eq. (54), one easily obtains 

dap _ n 2ap 1 + 4 + ^^/^ //3 + 4 ^^^^ •^Z' 

~dt ~ ~ ^^~c r^^^i ' 



-f- = -/34-(2e/3 + e/3sinV/3 + 2cos//3) 



dt c 6/3 



— = - ^Z2Z— (2 -e/3 cos //3) sin //3 , 



de/3 ^ Vm(i-/3)p/3 

dt r2 



(55) 



It is worth mentioning that dap /dt < for any time t. 

The set of differential equations Eqs. (55) has to be complemented with initial con- 
ditions. If the subscript "P" denotes orbital elements of the parent body (see Eq. (44)), 
then 

Pj-' 

Tp in = rp = - — ■ -f - ^PR 

1 + ep cos/p 

Vp in = Vp + A , (56) 
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where A is velocity with respect to the nucleus of the parent body, and 



_ P0 in 

in — ^ , 1. ^3 R in 

i- + e0 in COS //3 in 

in = R in ep Rin + Vp T in 6/3 



V/3 R in = \ e/S m Sm j„ 









in 


V (1 






in 



VpTin = \ (l + e/3i„COS//3i„) . 

V P/3 in 

ep Rin = (cos j„ COS Op in - sln f7;3 in sln „j COS ip in , 

sin flp in cosQpin + COS Qp in sinQp incosip in , sin 0/3 i„ sini^j„) , 

ep T in = (- cos ftp in sin Qp in - sin Qp in COS Qp in COS ip in , 

- sin Qp in sin 6/3 ,„ + cos ftp in cos Qp in COS ip in , COS 6/3 j„ sin ip in) , 

//3 in = Qfi in — in ■ (57) 

We write (epjv = epR x bpt) 

A = Aufl epR + Aut + At;/v cpjv • (58) 

Using Eqs. (44), (56) - (58), wc finally obtain (coniparc Gajdosfk and Klacka 1999): 
_ Pp {vt + Avrf + (Ati/v)^ 

P0 in - YZr^ ^2 ' 

f l + ep cos/p v^s _ y 

V 1-/3 4 y 

' 1 + ep cos/p UTs\^ / vr + At;fl \^ 
1-/3 ^^T / V '^T / ' 



p2 



Vt + . Aw/v ^ . . 

COS tp in = COSlp — COS6p Sinip, 

vts vts 



Vt + Avt . . , 

Simp in COS Up in = Sintp COSSip + 

Vts 



Av 



(cos6p cosip cosOp — sin6p sinfip) , 

Vt + Avt . . . p, , 
simp in sm Up in = sm^p sinSZp + 

Vts 

^ (cos6p cosip sin rip + sin6p cosfip) , 

VTS 

ep in COS fp in = (1 + ep cos/p) - 1 , 

Pp 

. _ 1 + ep cos fp Vts {vr + Avr) 

Gp in sm Jp in — _^ ^ ^2 5 

VT + Avt . . p, , 

sm I p in COS Qp in = Sinip COS6p + 

Vts 
Avn 

COSIp , 

Vts 

sin ip in sin Qp in = sinip sin6p , 

4s = {vt + AvtY + {AvNf , (59) 
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where vn and Vt are given by expressions presented in Eq. (78) ~ pp = ap (l 
For the special case A = Eq. (59) reduces to 



a,i„ = ap(l-/?)(l-2/?i±^£^) , (60) 



1 - el - 2(3 {I + ep cos fp) 
e,in = ^l ' (^1) 

where fp = Qp — u)p, ijjfj j„ has to be obtained from 

/3 + ep cos/p 



e/3 in cos (9p -U>I3 in) = 



l-f3 



6/3 in sm (Bp - W/3 in) = ^_ ^ , (62) 
^0 in = ) iff in = ip , 6/3 m = • (63) 

Physics of Eqs. (60) - (61) is following: meteoroid escapes from the Solar System when 
the orbital energy becomes positive and this can happen when the energy due to the 
radial component of the radiation force is included, without it being necessary for this 
force to exceed the gravitational attraction (Harwit 1963). Some figures may be found 
in Klacka (1992b). As an example we may mention that particle of /? = (1 — ep)/2 
moves in parabolic orbit if ejected at perihelion of the parent body, and, in an orbit with 
eccentricity = |1 — 3ep|/(l + ep) if released at aphelion; = for /? = ep = 1/3 
and aphelion ejection. 

By inserting Eqs. (55) into Eq. (50), taking into account that 6/3 i„ < 1 (see Eq. (61)), 
one can easily obtain for the secular changes of orbital elements: 

(^)=-/?^^^^. (64) 

4; 



dt' ^ c /. 2\'/' ' 

a/3 (^l-'-^-i 



«/3 

<f> = <^> = <^> = °- 

As an remark we may mention that a little more simple set of differential equations than 
that represented by Eqs. (64) (65) is obtained when semi-major axis ap is replaced by 
the quantity p/3 = ( 1 — e| ). It can be easily verified that a new set of differential 
equations for secular evolution is given by the following equations: 

i) dp^/dt = -2/3(/x/c)(l - 62)3/2/^/3, 
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ii) de^/dt = -(5/2)/3(/x/c)[(l - e2)3/2/p^]e0/p^. 

These two equations immediately yield p/s = p/sin {ep/epinY^^- 

The last relation corresponds to the relation between orbital angular momentum and 
eccentricity of the orbit: Hp = Hpin (e/j/e^m)^/^. 

Analogously, for secular evolution of osculating orbital energy (per unit mass) holds: 
Ep = vy2 - m(1 - P)/r ^ -m- el)/{2pp), 

Ep = Ep,4{l - e2)/(l - elJ] (ep.^jepfl^ 
Equations i) and ii) enable to write: 

i) dppldt = -2/3(m/c) [1 - e%Jpp/pp,nm'/' I PP. 

ii) dep/dt = -(5/2)/3(m/c) {eyi/plJ (1 - e^)3/2 / ef . 

These equations point out that near-circular / pseudo- circular orbits {ep « 0) and orbits 
with small values of pp are not described by the discussed secular changes of orbital 
elements and by the considered v/c approximation for the P-R effect. 

6.2. Secular changes of orbital elements - gravitation as a central acceleration 

As it was discussed in section 5, it is not wise to use — n {1 — /3i) e^j / as a central 
acceleration determining osculating orbital elements for the general case represented by 
Eqs. (40) or (41), since (3i changes almost randomly during a motion. In order to have a 
comparable results in disposal, we have to find secular changes of orbital elements when 
central central acceleration is given by — /x / r^. Thus, we will use — /x / as a 
central acceleration determining osculating orbital elements. 

Taking into account Eq. (49), we take the action of electromagnetic radiation as a 
perturbation to the two-body problem. We can immediately write for components of 
perturbation acceleration: 

Fh = /j4 - 2^4^, F, = -/3 4^, F, = 0, (68) 
r'' c r^ c 

where Fr, Ft and F/v are radial, transversal and normal components of perturbation 
acceleration, and two-body problem yields 

VRd = w - e sm / , 

V p 

VTd = (1 + ecos/) . (69) 

Perturbation equations of celestial mechanics yield for osculating orbital elements (a 
- semi-major axis; e - eccentricity; i — inclination (of the orbital plane to the reference 
frame); f2 - longitude of the ascending node; u - longitude of pericenter; Q - Q = oj + f 
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is the position angle of the particle on the orbit, when measured from the ascending node 
in the direction of the particle's motion): 



da 
rfi" 
de 
di 
di 
di 

dn 

m 

dw 



2 a 



1 - e2 

fp 



^ {Frc sin f + Ft {1 + e cos/)} 



^ {Fr sin f + Ft 



Fjv cos 6 
sinG 



cos/ 



cos/ 



1 + ecos/ 



/MP 
r 

/JPp 
1 



Fn 



smi 



IP 



e 
r 



dQ 
dt 



sin e 

tn —. — r COS I 



2 + e COS / 
1 + e COS / 



sin/ 



smi 

r 



sine 

-TAT ^ — r COS J 



V P ^^-"^ * 
where r = p/{l + e cos /). 

Inserting Eqs. (68) - (69) into Eq. (70), one easily obtains 

da 



(70) 



2/3 



dt 
de 
di 
di 
di 

dn 

'dt 
duj 
'dt 
dB 
'dt 



1 



. . ^ u 2 a 1 + + 2e cos / + sin^ / 
2 e sm/ - '-^ 

z c 1 — e 



^ sm/ 



/? (2e + esin2/ + 2cos/) 



= 0, 
= 0, 



/3 



fjpp 1 



cos/ - /3 4 (2 - ecos/) sin/ , 



/HP 



(71) 



It is worth mentioning that da/ dt < for any time t docs not hold (perturbation corre- 
sponds to complete nongravitational acceleration, and, thus, Eqs. (55) do not hold). 

The set of differential equations Eqs. (71) has to be complemented with initial condi- 
tions. If the subscript "P" denotes orbital elements of the parent body, then the initial 
orbital elements for the particle ejected with velocity A are: 
{vt + Avt)"^ + (At;jv)^ 



Pi: 



Pp 



(1 + ep cos/p) 



cos li 



sini,;„, cos flir, 



(1 

Vt 



ep cos/p) 



2 

vts] ^ / Vr + Avr \ ^ 



Vts 
Vt + A( 



■ cos Ip 



Av 



Vts 



Vts 

— sin i p cos Hp - 



cos9p sinzp , 
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— (cosQp cosip cosflp — sinOp sinfip) , 

VTS 

■ n '"T + AvT . . . o , 
smiin smilin = simp smilp + 

VtS 

^ (cosOp cosip sinJlp + sin Op cos ftp) , 

VtS 

Bin COS fin = — (1 + Cp COS/p) - 1 , 
PP 

Bin Sm fin = (1 + ep COS/p) 2 ' 

r, + . . ^ , 

sm^j„ cosBi„ = sm^p cos Bp + 

COSlp , 

VTS 

siniin sin6i„ = sinip sinOp , 

4s = {vT + Avrf + (AvNf , (72) 

where vr and vt are given by expressions presented in Eq. (44) - pp = ap (l — ep). 

For the special case A = Eq. (72) reduces to a simple fact: initial osculating orbital 
elements of the particle are identical with those of the parent body: 

^in — ^P ) ^in — ^P ; ^in — ^P j 

^in = ^P ) l^in = t^P j ©m = • (73) 

The important fact is that Eqs. (71) contain also terms not proportional to v/c (<C 1). 
These important terms protect us to use procedure analogous to that represented by Eq. 
(50). While dispersion of osculating orbital elements is very small during a time interval 
T for the case when fj, {1 — (3) is used in central acceleration, the dispersion of osculating 
orbital elements may be large during the same time interval for the case when fj, is used in 
central acceleration (compare Figs. 1 and 2 in Klacka 1994b). Thus, any formal averaging 
of Eqs. (71) leading to equations analogous to Eqs. (64)-(67) is not correct. 

We have explained that it is not allowed to make a simple time averaging analogous to 
that described by Eq. (50), when — fj, br / r"^ is used as a central acceleration determining 
osculating orbital elements. However, it is of interest if we have to numerically solve 
Newtonian vectorial equation of motion (Eq. (49)) and make numerical time averaging 
(over a time interval between passages through two following pcriccntcrs) , or if some 
analytical simplifications can be done something analogous to Eqs. (64)- (67). 

Fortunately, wc can make analytical calculations for the purpose of obtaining secular 
changes of semi-major axis and eccentricity even when — /i ep / is used as a central 
acceleration determining osculating orbital elements. We will derive the correct equations 
in the following two subsections. 
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6.2.1. Radial forces and orbital elements 

We will proceed according to Klacka (1994), in this subsection. 
Let us consider a gravitational system of two bodies 

v = -^en. (74) 

Let perturbation acceleration exists in the form 

« = /3 ^ efl , (75) 
< /3 < 1. Thus, the final equation of motion is 

* = -^^-. ,™, 

Eq. (76) yields as a solution the well-known Keplerian motion and the orbit is given by 

(77) 



1 + Be cos (6 — UJc) ' 

where 



p, = ae (1 - el) . (78) 

The subscript "c" denotes that orbital elements are constants of motion. If we write 

V = VcR en + VcT ex , (79) 

where bt is a unit vector transverse to the radial vector eR in the plane of the trajectory 
(positive in the direction of motion), we have 



VcR = Vf^ (1 - P) Pc ^ ec sin(e - Wc) , (80) 



VcT = (1 - /3) pc^ [1 + Cc cos (6 - Wc)] . (81) 

In principle, we may consider also a new set of orbital elements, which are defined by 
the central gravitational acceleration. Eqs. (77) - (81) are then of the form 

^ (82) 



1 + e cos (6 — u>) ' 

p = a (1 - e^) , (83) 

V = VdR Br + VdT er , (84) 

VdR = a/m e sin (6 - w) , (85) 



VdT = [1 + e cos (G - uj)] ; (86) 
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the fact that is unchanged in both sets of orbital elements is used; br = ( cos 9, sin© 
), bt = { — sinG,cos0 ). 

Position vector and velocity vector define a state of the body at any time. Equations 
(77) and (82) yield then 

^ P _ (87) 

1 + Cc cos (0 — ojc) 1 + e cos (9 — w) ' 

Analogously, the other two pairs of equations (Eqs. (80) and (85), and, Eqs. (81) and 
(86)) give 

^J{1 - (3) pc^ Cc sin (9 - ujc) = \^ e sin (9 - w) , (88) 



y^(l - /?) pc^ [1 + Cc cos (9 - LOc)] = \fp^ [1 + e cos (9 - w)] . (89) 
One can easily obtain, using Eqs. (87) and (89), 

Pe (1 - /3)=P, (90) 
and Eqs. (88)-(89) yield then 

(1 - /3) Cc sin (0 - Wc) = e sin (9 - w) , (91) 
(1 - /3) [1 + Cc cos (9 - Wc)] = 1 + e cos(9 - w) . (92) 
Eq. (92) may be written in the form 

(1 - /3) ec cos (9 - Wc) - /3 = e cos (9 - w) . (93) 
Eqs. (91) and (93) yield 

= (1 - fif el + (3^ - 2 p {1 - P) ec cos (9 - Wc) • (94) 

Equation for w is given by Eqs. (91) and (93), using also Eq. (94). Finally, Eqs. (78), 
(83), (79) and (94) yield 



/i J. ^ 1 + eg + 2ee cos(9 - Wc) \ ' 
a = ac |1 + P 1 - el j ' ^ ' 

Eqs. (94)-(95) show that orbital osculating elements a and e change in time, during an 
orbital revolution - the larger /3, the larger change of a and e. 

The osculating orbital elements e and a obtain values between their maxima and 
minima, which can be easily found from Eqs. (94)-(95). One can easily verify that these 
relations hold: 

emin = I (1 - /3) Cc - /3| < e < (1 - /?) Cc + /3 = Cmax , (96) 
ttc 1 - Be + /3 (1 + Cc) ~ ttc ~ 1 + Cc + /3 (1 - Cc) ttc 
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6.2.2. Mean values of semi- major axis and eccentricity 

Eqs. (96) and (97) represent interval of values for eccentricity and semi-major axis, when 
— fj, cr / r'^ is used as a central acceleration determining osculating orbital elements. 
However, we can make time averaging during a period T, which was decribed by Eq. 
(50): 

(e) = ^ ^ e{t) dt, {a) = ^ a{t) dt , 

r VM U P)Pc, ' ^ 1 + ee cos/e • ^ ' 

Eqs. (98) yield 



e(/c) 



(1 + ec cos/c) 

2n 



dfc , 



^ TT ^0 (1 + Cc COs/c) 

Using Eqs. (94) and (95), we finally obtain 



(99) 



, , 2,3/2 1 /-^^ ^(l-/?)'e2+/32_2/3(l-/3)eeCOS/e 

(e) = (1 - ej — / — —2 dfc , 

^0 (1 + ec COS fc 



(l-t-Cc cos/c 



I -1 



<„, = eS='^ ± [1. (1 + eg + )/(!-.;)] ^^^^^ 

^0 (l-hCc cos/c) 

The following properties can be verified: 
i) (e) > P, (e) = p for = 0; ii) (e)/ec > 1, (e) = for /? = 0; 
iii) d{e)/dec > 0; iv) d{e)/df3 > 0; 

v) (a) > ac/(l + /3), (a) = ac/(l + /?) for Cc = 0; 

vi) {a) /(ic < 1, (a) = flc for /3 = 0; 

vii) d{a)/dea > 0; viii) d{a)/dac > 0; ix) d{a)/dp < 0. 

6.2.3. Secular changes of semi-major axis and eccentricity 

Summarizing our results, it is possible to calculate secular evolution of eccentricity and 
semi-major axis, according to the following prescription. 
At first, initial conditions for ap and are calculated: 



M.„ = -p(i-^)(i-2g 'V-T'' ) 



supposing that particle was ejected with zero ejection velocity from a parent body - 
quantities with subscript "P" belongs to the parent body's trajectory; as for more general 
case, Eqs. (56), (58) and (59) have to be used, {a^)^ = P0 in / (1 - j„), (e/3)m = ^f} in- 
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As the second step, the set of the following differential equations must be solved for 
the above presented initial conditions: 
2 + 3e? 



rrz = _ /3 ^ /3 



dt c /, ' 



Finally, semi-major axis and eccentricity are calculated from: 



a 



< 2^3/2 1 Z-^- [l + /?(l + e^ + 2e;3Cosar)/(l-e^) 
^0 (1 + Cflcosa;) 



(1 + cosa;)^ 
3/2 1 r^- ^(l-/3)%2+/32_2/3(l-/J)e0cosa; 



e = (1 - e|)^/^ ^ / V^-^^ e,.p -.p,.-p;e,cos. _ 

^ 27r io (l + e/3Cosa;)^ ^ ^ 

The set of equations represented by Eqs. (101)-(103) fully corresponds to detailed 
numerical calculations of vectorial equation of motion, if we are interested in secular 
evolution of eccentricity and semi-major axis (supposing (e/3)i„ < 1 and e/j does not 
correspond to pseudo- circular orbit) for the case when central acceleration is defined by 
gravity alone. 

It is worth mentioning that instantaneous time derivatives of semi-major axis and 
eccentricity may be both positive and negative, while secular evolution yields that semi- 
major axis and eccentricity arc decreasing functions of time. 

We may mention that a little more simple procedure is obtained when semi-major 
axis a/3 is replaced by the quantity = ( 1 — e| ). As it was already mentioned 
behind Eq. (67), 

i) dv^ldt = -2/?(m/c) [1 - ej ,Jpp/p0 ,„)5/2]3/2 / p^, 

ii) dep/dt = -(5/2)/?(m/c) {e'Jl/p} ,J (1 - e^)^/^ / ef . 

Initial condition for pjs is given by the first relation of Eq. (59), e. g., Pjs m = Pp / (1 — /?) 
for zero ejection velocity. The first integral of Eq. (103) is replaced by a simple relation: 
p = P0 {1 - d) (Eq. (90)). 

(It is important to stress that quantities p = {p), and a = (a), e = (e), present in Eq. 
(137), do not fulfil relation p = a{l- e^).) 

7. P-R effect — equation of motion to the second order in v/c 

We can write, on the basis of Eqs. (48) (Balck and Klacka 2002) 



dt r'' 
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where fi = G Mq, cr = e\ and /? = /3iis a non-dimensional parameter ("the ratio of 
radiation pressure force to the gravitational force"; see also Eq. (43)). The first term in 

Eq. (104) is Newtonian gravity for two-body problem, the second term is the Poynting- 
Robcrtson effect in flat spacctimc, the third term corresponds to Einstein's correction to 
Newtonian gravity and the last term is a sort of interference term describing the mixing 
between the effects of gravity and radiation pressure. 

7.1. Secular changes of orbital elements - radiation pressure as a part of central acceler- 
ation 



We have to use — /U (1 — /?) e/j / as a central acceleration determining osculating 
orbital elements if we want to use a fact that the elements do not change rapidly during 
a motion described by Eq. (104) - during a time interval T, where T is time interval 
between passages through two following pericenters / perihelia: 

dv /i(l — /3) _ ^ fi f\v-eii 

dt r' 
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We can immediately write for components of perturbation acceleration to Keplerian mo- 
tion, on the basis of Eq. (105) - /3 is considered to be a constant during the motion: 
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(106) 



where r, Fjs t and ^ are radial, transversal and normal components of perturbation 
acceleration, and the two-body problem yields 



Vf3 R = 
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//x(l- 
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P/3 



e/3 sin fi3 , 

{l + e^cosf^) , 



(107) 



where ep is osculating eccentricity of the orbit, fp is true anomaly, pp = ap{l — e^), ap 
is semi-major axis. The important fact that perturbation acceleration is proportional to 
v/c (<C 1) ensures the above metioned small changes of orbital elements during the time 
interval T. 
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Perturbation equations of celestial mechanics yield for osculating orbital elements (a/3 
- semi-major axis; Cfs - eccentricity; ijs - inclination (of the orbital plane to the reference 

frame); Qp longitude of the ascending node; LUp longitude of periccntcr / perihelion; 
0^ is the position angle of the particle on the orbit, when measured from the ascending 
node in the direction of the particle's motion, Qf3 = ojp + fp): 
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where r = pp/{l + ep cos fp). 

Inserting Eqs. (106) - (107) into Eq. (108), one easily obtains 

,2 



dap 
~dF 



_ /X 2a/3 1 + + 2ep cos fp + e| sin^ fp 

~ P~J} I 19 "I" 



M ap Ui{l-P) pxXai + 6xX, 



a2 



1 - 



14 
1~ 



Pf3 
ep sin fp 



1 - p2 

1 Bp 



1 - 



1-/3 



isin(2//3) + ep sin fp 



Xa2 

dep 
dt 

dip 
dt 

d£lp 

"df 
diOp_ 

dt 



(1 + r=^) //5 + (1 + r=^) 4 sin (2//3) + e| sin fp 



13^ ^ {2ep + ep sin^ fp + 2 cos fp) + 



fi 1 Ui{i-p) (I 



Pl3 



14 
4 



sin//3 + 1 



sin//3 



7 



1-/3 
2 



1 + 

0, 
0, 



(2 -6/3 cos //3) sin //3 

cep 



ep sin (2//3) + 36^ sin//3 
6/3 sin (2//3) + 7e| sm fp 



(108) 



29 



6/3 y P/3 \2 
Xu,2 = T^^) ^^'^ - ^6 + e/3 cos^ //3 + e| cos fp 



(109) 



We want to find secular changes of osculating orbital elements up to 
As for the terms proportional to 1/c^ in Eq. (109), we may take a time average (T is 
time interval between passages through two following pericenters) in an analytical way 

Vf^ (1-/3) 1 r i.^fdfp^-' 



\/m (1 - /3) 1 r2 



■/o ^^^^^ 



dU 



^3/2 27r x/m (1 - /?) 



assuming non-pseudo-circular orbits and the fact that orbital elements exhibit only small 
changes during the time interval T; the second and the third Kepler's laws were used: 
dfjj/dt = -\//I(r^^"^)p^ - conservation of angular momentum, a^^/T"^ = /x(l — /3)/(47r^). 
The result is: 

^ItF^'' - ° ' 
^^^^^ - °' 

- c^af (l-e^) (1^ • ^'"^ 

Taking into account terms proportional to 1/c in Eq. (109), it is inevitable to use 
perturbation theory of the second order. Simple averaging represented by Eq. (110) is 
not sufficient: the last equation of Eq. (109) yields 
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Eq. (146) holds due to the fact, that Fp jv = - in reality, equation for dQfj/dt presented 
in Eq. (108) has to be used. Instead of Eq. (110), we have more precise method of 
averaging, now. We can write, using Eq. (112): 
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We are interested in secular changes of orbital elements up to the order l/c^. Within 
this accuracy we can write for the terms on the right-hand side of Eq. (113) 
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We have to find osculating orbital elements present on the right-hand sides of Eq. 
(109) in terms of true anomaly fp to the required accuracy in l/c^. Differentiation of the 
relation = { 1 — ) yields, using Eq. (109) for dap/dt and de/}/dt, 
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and, using Eqs. (109), (112) and (114), we obtain 
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Now, if we restrict ourselves to the first order in 1/c, Eqs. (117)-(118) reduce to 
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Eqs. (119)-(120) can be easily solved: 



e/3 = 6/30 - 



2sin//3 + 
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5//3 - - sin 2//3 
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(122) 



where initial values PfM, e^o correspond to fp^ = 0. 

On the basis of Eq. (121), we can immediately write solution of Eq. (117): 
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On the basis of Eqs. (121) - (122), we can rewrite Eq. (118) to the form 
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where 6/30 has to be inserted into expressions for X^i and (see Eq. (109)) instead of 
6/3. Solution of Eq. (124) is: 
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(125) 



Semi-major axis a/3 can be obtained from the relation p/3 = a/3 ( 1 — e| ), using Eqs. (123) 
and (125). 
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Let us calculate secular change of u)0 to the order generated by osculating orbital 
change proportional to 1/c ~ compare Eqs. (109) and (HI). Putting Eq. (114) into Eq. 
(113), we can write to the required accuracy, 



since time derivative of is proportional to 1/c, and, thus, the third term on the right- 
hand side of Eq. (148) would yield a higher order than in the last integral. 
The last integral in Eq. (126) consists of two terms: 
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Inserting relation for dup/dt from Eq. (109), we can immediately write for I^i, within 
the required accuracy: 
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where results of Eqs. (121) and (122) were used, and, Hfj = ^ — f3)pp. 
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As for calculation of /ai2i, we need to consider that Hfj = sj /i(l — /3)p/3 and Eq. (121). 
We obtain: 
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Calculation of Ioj22 is analogous to calculation of Iui2i, but we have to use Eqs. (121) and 
(122) simultaneously. The result is: 



13"^ H Att ( 7 1 



'^^^ 1 - /? c2 pf^o \ 2 e/30 
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Eqs. (129)-(131) yield 



1- (3 POQ y 8 2 e/30 y 
Eqs. (127), (128) and (132) yield 
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We have calculated the last integral in Eq. (126) - it is represented by Eq. (133). 
Since it is proportional to 1/c^, the first integral in Eq. (126) has to be proportional to 
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Finally, Eqs. (126), (133) and (134) lead to 

(^)r - ^ (-'- + ^] (135) 

Total shift of pcricentcr / perihelion is given as a sum of Eqs. (Ill) and (135). Chang- 
ing the index /30 into /3 and using p^ = (1 — e|), we finally receive: 
d^ _ 3/x3/2 i + /?2 (_i3/8 + 7/e^)/3 

It can be easily verified that < dLOf} / dt > is 

i) an increasing function of /?, 

ii) the perihelion circulates in a positive direction, 

iii) the rate of the advancement of perihelion is not bounded for /3 — > 1. 

Let us calculate secular change of ap to the order 1/c^, generated by osculating orbital 
change proportional to 1/c compare Eqs. (109) and (111). Putting Eq. (114) into Eq. 
(113), we can write to the required accuracy, 
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since time derivatives of ap and ivp are proportional to 1/c, and, thus, the third term on 
the right-hand side of Eq. (114) would yield a higher order than in the last integral. 

The last integral in Eq. (137) consists of two terms: 
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Inserting relations for dap/dt and dwp/dt from Eq. (109), we can immediately write for 
lai, within the required accuracy: 
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where results of Eqs. (121) and (122) were used. 

2ap (l + el + 2ep cos fp + sin^ fp 



'o2 



- r 

c Jo 



H, 



(l - el) 



dfp. 



(140) 



Considering that Hp = a/^(1 — P)pp and Eqs. (121) and (122), we have 
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Inserting results of Eq. (141) into Eq. (140), one finally obtains 
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TT + (IOtt + 3) e|o + (Qtt + 2) e^Q ' 



(142) 



We have calculated the last integral in Eq. (137) - it is represented by Eq. (142). 
Since its dominant part is proportional to 1/c, the first integral in Eq. (137) has to be 
proportional to 1/c^. Thus, the first integral in Eq. (137) is given, according to Eq. (115), 
as 
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It can be easily verified that the last integral in Eq. (143) equals zero, within the required 
accuracy, and 
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Finally, Eqs. (137) - (140), (142) - (144) yield 
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Total secular change of semi-major axis is given as a sum of Eqs. (Ill) and (145). 

Changing the index /30 into (3 and using Eq. (144), we can finally write 
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(1 + 6/3 cos x) 27r Jo (1 + 6/3 cos a;) 

Let us calculate secular change oiep to the order l/c^, generated by osculating orbital 
change proportional to 1/c compare Eqs. (109) and (111). Putting Eq. (114) into Eq. 
(113), we can write to the required accuracy. 



36 



since time derivatives of e^? and LOf^ are proportional to 1/c, and, thus, the third term on 
the right-hand side of Eq. (114) would yield a higher order than in the last integral. 

The last integral in Eq. (147) consists of two terms: 

- io H^Hf, dt dt ""^^^ 
_ de, 
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Inserting relations for dep/dt and dujjs/dt from Eq. (109), we can immediately write for 
/el, within the required accuracy: 
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where results of Eqs. (121) and (122) were used. 
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Considering Eq. (122) and the first relation of Eq. (141) for I/H/3, one finally obtains 
from Eq. (150): 
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Putting the results represented by Eqs. (144), (148), (149) and (151) into Eq. (147), 
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Total secular change of eccentricity is given as a sum of Eqs. (Ill) and (152). Changing 
the index (30 into /3 in Eq. (152), we can finally write 
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Let us calculate secular change of pp (see Eq. (116)) to the order 1/c^, generated by 
osculating orbital change proportional to 1/c - one can easily verify that (dpp/dt)ii — 0. 
Taking into account Eqs. (109) and (116), wc can write to the required accuracy, 
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since time derivatives of pp and ujp arc proportional to 1/c, and, thus, the third term on 
the right-hand side of Eq. (114) would yield a higher order than 1/c^ in the last integral. 
The last integral in Eq. (154) consists of two terms: 
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Inserting relations for dpp/dt and dujs/dt from Eqs. (116) and (109), we can immediately 
write for Ipi, within the required accuracy: 
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where Eq. (121) was used. Putting results of Eqs. (143), (144), (155) and (156) into Eq. 
(154), one obtains 
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Total secular change ofp/j is given as a sum of Eqs. (157) and {dpi3/dt)ii = 0. Changing 
the index /30 into (3, wc can finally write 
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Comparing Eqs. (153) and (158), one easily obtains 
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This can be easily integrated and the results may be written as: 
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8. Secular change of advancement of pericenter /perihelion — gravitation as a 
central acceleration 

We will use — /U e/j / as a central acceleration determining osculating orbital elements. 
As we have seen in the previous subsection, corrections of the order {v/cf' represent very 
small corrections with respect to the order v/c - only dujjs/dt is important, as for secular 
changes. Thus, the most important results, for the case when — fj, en / r"^ is used as 
a central acceleration, were presented in section 6.2. However, we have not calculated 
secular change of dco/dt in the section 6.2. We will do this now, since results of the 
preceding section 7.1 have to be used. 

We will use Eqs. (91), (93) and (94) in the form 
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We admit that the subscript "c" may be changed into /?. Using Eq. (162), one immediately 
obtains 
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Inserting the right-hand side of the second of equations of Eq. (161) into Eq. (163), one 
finally, after differentiation, obtains 
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Differentiation of the third equation of Eq. (161) yields 
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Inserting Eq. (165) into Eq. (164), we can write 
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Eq. (166) is the decisive equation which enables us to find secular change of w. 
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As a first approximation, let us consider Keplerian orbit characterized by conditions 

dttc/dt = dpc/dt = dcc/dt = diOc/dt = 0, Pc = ctc(l ^ e^). Using averaging of the type of 
Eq. (98), or, Eq. (110), we can write for Eq. (166) (we use r'^dfc/dt = ^/]I{r^^'p)pc) 



(-) - - / 



dt ' 



dco 
Itt 



dt 



2-n- 



Jo 



doj , 2 ij- 1 1 



(1 - f])' el + 13^- 2/3 (1 - (3) e„ cos /, 



v/m(1-/?)Pc d/c • (167) 



If the result ^ dx/{l + A: cos a;) = 27r/-\/l — k"^ is used, one finally obtains 
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if/3<ee/(l + ee) 

^IJ,{l-P)/al'^/2 if /3 = ee/(l + e,) ■ 



(168) 



3/2 



if /3 > e,/(l + ee) 



Let us remind that dOc/rft = n{l - (3) / aj"^ (see Eq. (100)). 

Let us consider P-R effect to the first order mv/c. We will use Eq. (55) and also Eq. 
(112), which we summarize in the following Eq. (169): 



de0 
dt 
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where = -^/^(l — P)p0. Changing subscript "c" into "/3" in Eq. (166) and inserting 
the third of Eq. (169) into Eq. (166), we can write 



doj 
'dt 



1 + 



/3'-(l-/3) e 



2 ^2 



/3'-(l-/3) e 



2^2 



dcof, 
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= (1 - pf 6l + - 2p {1 - /3) 60 cosf0 . 



(170) 



On the basis of Eqs. (113)-(115), (143)-(144), (169)-(170), we can write within the re- 
quired accuracy: 

r2ir ^2 / ^2 



.dev. , r2 / r2 du0\ dw 

= ^ / 
2 io 



1 + 



(1 - /?)' el + p^- 2/3 (1 - /3) 60 cos //, 
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27ra' 



730 

2ir 
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(171) 
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Using also Eq. (122) in the integral of the second line of Eq. (171), we finally receive 



^dt^ ' 3/2 \ 



a/3 



^ = 



(1 - e$f' [5 73 (e/5) - 2 h (e^)] (p - 




_/32 + (l-/3)2e| 

2/3(1-/3) 
(l-/3)%|+/32 ' 



1 X 



Ia{e) ^ TT 77- v?'^^;, a = l,2,3, (172) 

27r Jo (l + ecosx) 

where ■i?ij(a;) = 1 if a; > 0, t?ij(a;) = if a; < (Heaviside's step function); it is assumed 
that /3 ^ 6/3/(1 + 6/3). We see that advancement of pericenter /perihelion exists already in 
the first order in u/c if gravity alone is taken as a central acceleration. 

9. P-R effect and near circular orbits 

When the orbits arc near circular (pseudo-circular) - central acceleration contains ra- 
diation pressure term the orbit can not reduce in semi-major axis without increasing 
in eccentricity. Both types of orbital elements, defined by central accelerations, have 
been used in detail numerical calculations in papers by Klacka and Kaufmannova (1992, 
1993). Due to the property (e) > (3 (see section 6.2.2), one must be aware that also values 
(e) > /3 may correspond to pseudo-circular orbits. The results for pseudo-circular orbits 
were analytically confirmed by Breiter and Jackson (1998). As an advantage, the analyt- 
ical approach reproduces known results without detail numerical calculations. However, 
the analytical approach paralelly produces nonphysical results which may not be distin- 
guishable from the correct results. The nonphysical analjdiical results are caused by use 
of the P-R effect in the first order in u/c - very special analytical solutions will diminish 
when higher orders in u/c are used. The nonphysical analytical results have been discussed 
in more detail by Klacka (2001; see also http://xxx.lanl.gov/abs/astro-ph/0004181). 

10. Solar wind effect 

We have calculated secular changes of orbital elements for the P-R effect up to the second 
order in v/c. The effects of the second order in w/c seem to be small to play an important 
role in Solar System studies. In reality, similar effect coming from the Sun exists and 
this effect may play more important role. Although this effect is not connected with 
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electromagnetic radiation, its secular changes of orbital elements to the first order in 
v/c (more correctly, v/u, where u is the speed of solar wind particles) correspond to 

the secular changes for P-R effect. This effect is caused by solar wind particles hitting 
an interplanetary dust particle. The aim of this section is to obtain secular changes of 
orbital elements of the interplanetary dust particle under the action of the solar wind, 
up to the second order in v/u. 

Let us consider equation of motion in the form (we neglect decrease of particle's mass) 



d V. 



dt 



l3 fj, u 



{[(> 



V ■ Sr 



u 



eR 



+ 



(v-eR\ V 1 r/v\2 /v-ei?\2l 1 



(173) 



where, as standardly in this paper, ii = G Mq and non-dimensional parameter /3 is 
defined in Eq. (43) ("the ratio of radiation pressure force to the gravitational force"); 
rj « 1/3. Adding the right-hand side of Eq. (49) to the right-hand side of Eq. (173), one 
obtains 

dv _ jj, ^ ^^^/('x ^'^■'^^e ^]- + 
d t r^lV c J cJ 



Q[ r2 c llv u J 
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V 1 
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eR 
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(174) 



The first term in Eq. (174) is Newtonian gravity for two-body problem, the second term is 
the Poynting- Robertson efii'ect in flat spacetime, the third term corresponds to the action 
of the solar wind up to the second order in v/u. 
Eq. (174) may be written in the following form: 

dv Lb ^/ n u\ Li ^ / rj\ a fv-BR v\ 
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c u 



{V ■ BRf 



c u 



(175) 



10.1. Secular changes of orbital elements - radiation pressure as a part of central accel- 
eration 

Neglecting solar wind pressure term ri{fi/Q[){fi/r'^){u/c), we can rewrite Eq. (175) to the 
following form: 

dv _ /u(l-/?) 



dt 



eR 



/34 



eR , V- 

C Cj 
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Bt 



(176) 
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where the second term causes deceleration of the particle's motion. We use 
— ^ (1 — /3) / as a central acceleration determining osculating orbital elements; 

P is considered to be a constant during the motion. On the basis of Eq. (176), we can 
immediately write for components of disturbing acceleration to Keplerian motion 



Q'l J c Q[ r'^ cu ' 



,2 



jp _ r.l.,V\lJ'Vf3T filJ-]vf)RVpT,l^l3T 



Q[Jr'^c Q'l ''^ 1 CM ' 2 cu 

F0N = 0, (177) 

where Fp ji, t and Fp n arc radial, transversal and normal components of perturbation 
acceleration, and the two-body problem yields 



Vf3 R = \ e/3 sm ffj , 



If, (1- 


/3) 




//i (1 - 


0) 



V0T = \ — ^ (l + e/3Cos//3) , (178) 

V P/3 

where Cjs is osculating eccentricity of the orbit, is true anomaly, Pjs = a/}{l — e|), Ujs 
is semi-major axis. The important fact that perturbation acceleration is proportional to 
v/c ("C 1) ensures small changes of orbital elements during a time interval T [T is time 
interval between passages through two following pericenters / perihelia). Thus, we may 
take a time average in an analytical way 

(..^^r.<o-^ir./..(f)-'* 



(1 - P) 1 
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/■27r ^2 

Ja Ja{l-B)pB 
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9U0) df0 , (179) 



2 27r 



assuming non-pseudo-circular orbits and the fact that orbital elements exhibit only small 
changes during the time interval T; the second and the third Kepler's laws were used: 
dffj/dt = \/ iJi{l — P)p/3 - conservation of angular momentum, a^/T'^ = — 
Relevant perturbation equations of celestial mechanics yield for osculating orbital 
elements (w^ - longitude of pericenter / perihelion; 6/3 is the position angle of the particle 
on the orbit, when measured from the ascending node in the direction of the particle's 
motion, 6/3 = ujp + Jb)- 



dt 1 - e| Y /i (1 - /3) 
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dujf-. 



1 
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17 f 2 + (' J cos jj3 . . 

1 + 6^3 cos ' 



(180) 



dt ep V m(1 - 

where r = Pf}/{1 + cos f^). 

Putting Eq. (177) and (178) into Eq. (180), making procedure of averaging of the type 
given by Eq. (179), we finally receive 
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Initial conditions are given by Eq. (59), or, Eqs. (60) - (62). 

Similarly, for the quantity p/3 = ^1 — e^^ , one can easily obtain 
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(181) 



(182) 



We come to the conclusion that the following relation results from Eqs. (181) and 
(182): 



-4/5 -4/5 
P/3 ' = P/3 m /„ , 



(183) 



as for secular changes of and e/3, for the simple case of equation of motion of an 
interplanetary dust particle under the action of solar wind represented by Eq. (173) and 
for the P-R effect. Thus, Eq. (182) can be solved as a separate equation 
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10.2. Secular changes of orbital elements - gravitation as a central acceleration 



We can immediately write, on the basis of Eqs. (103), (172) and (181): 
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a = a^(l-e^) -^^ 
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and, for secular change of longitude of pericenter/perihelion we have 
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(188) 



(1 + ecosx) 

where 'dnix) = 1 if x > 0, ^h{x) = if a; < (Heaviside's step function); it is assumed 
that P ^ 6^/(1 — ep). 

Initial conditions arc given by: 

i) Eq. (59), or, Eqs. (60) (62) for the set of Eqs. (186) - (187), and, 

ii) Eq. (59), or, Eqs. (60) - (62) and Eq. (72), or, Eq. (73) for the set of Eqs. (186) and 
(188): Eqs. (72) - (73) are required for initial value of oj. 



10.3. Solar wind - discussion 

The results presented in this section hold for the most simple approximation of the solar 
wind action - only radial component of the solar wind particles is considered (Eq. (173) 
is taken as an approximation to more general equation of motion presented in Klacka 
and Saniga 1993; see also Leinert and Grun 1990). Moreover, solar wind causes decrease 
of particle's mass and the secular change of particle's mass m (present in /3, see Eq. (43)) 



is given as dm/dt = —KA'^^j / (a^y 1 — 6^), where K is n constant depending on the 
material properties of the particle and ^'.jy is the proper effective cross sectional area of 
the particle. If A'^^ ^ (area toward the Sun) is changing during the particle's motion, one 
has to use dm/dt = —KA'^j^il — vp jilu^jr^ and no averaging is possible considerations 
made in section 16 hold for the case ^g^^ = where A! was defined above Eq. (16), and, 
moreover, ^gjj does not change during the particle's motion (e. g., spherical particle). 
Real velocity vector of solar wind particles is nonradial and the nonradial component 
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increases with decreasing distance from the Sun (e. g., Stix 2002). As a consequence, real 
solar wind effect may cause acceleration of meteoroids in small distances from the Sun, 

instead of their deceleration. 

Finally, solar wind causes also charging of meteoroids and Lorentz force has to be 
taken into account in the case of submicron grains. 

11. Summary and conclusions 

The paper derives and presents relativistically covariant equation of motion for dust par- 
ticle under the action of electromagnetic radiation - see Eq. (40). As for most frequent 
applications to systems in the universe (e. g., meteoroids in the Solar System, dust par- 
ticles in circumstellar disks), equation of motion in the form of Eqs. (41) and (42) are 
sufficient: application of Eq. (41) (for the case F^j = 0, j = 1, 2, 3 and under some 
assumption about particle's rotation) may be found in Kocifaj et al. (2000), Klacka and 
Kocifaj (2001). Some other accelerations may be added to the right-hand side of Eq. (41) 
~ e. g., gravitational perturbations of planets, solar wind effect (sec Eq. (173) or some 
more precise form of equation of motion) or some other nongravitational accelerations. 

Special attention was devoted to the Poynting-Robcrtson effect, since this effect is 
standardly used in Solar System studies. We have derived secular orbital evolution for 
the P-R effect up to the second order in v/c. General equation of motion for interac- 
tion between particle and incident electromagnetic radiation shows that radiation cannot 
be considered as a part of central acceleration. The central acceleration has to contain 
only gravity of the central body (star/Sun); moreover, this corresponds to the physical 
situation when radiation effect is considered as a disturbing effect. In order to compare 
secular changes of semi-major axis and eccentricity for real cosmic dust particle and the 
P-R effect, the paper derives and presents also secular changes of these orbital elements 
for the P-R effect: see Eq. (103) in section 6.2.3 and Eq. (172) in section 8 - advancement 
of pericenter/perihelion exists even in the first order of v/c. 

Solar wind effect is also considered in section 10. Solar wind, in its simplest approxi- 
mation, produces secular changes of the orbital elements analogous to the P-R effect, as 
for the first order in v/c see first two equations in Eq. (181) and Eq. (182). However, 
second order in solar wind effect produces more significant changes of orbits than it is in 
the case of the P-R effect, if radiation pressure is a part of central acceleration compare 
Eqs. (136), (146), (153) and (181). Section 10.2 presents secular changes of orbital ele- 
ments for the P-R effect and the most simple approximation of the solar wind effect when 
solar gravitation alone is considered to be a central acceleration, in its simplest produces 
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Application to larger bodies, e. g., asteroids, may be found in Klacka (2000c) - some 
kind of thermal emission has to be added (quantities Q'^j, j = 1, 2, 3 present in Eqs. (40) 
- (42) have to be calculated). 

Appendix A: Another formulation of the equation of motion 

(Reference to equation of number (j) of this appendix is denoted as Eq. (A j). Reference 
to equation of number (i) of the main text is denoted as Eq. (i).) 

Proper reference frame of the particle - stationary particle 

The equation of motion of the particle in its proper frame of reference is taken in the 
form 



where E' is particle's energy, p' its momentum, r is proper time, S' is the flux density of 
radiation energy (energy flow through unit area perpendicular to the ray per unit time), 
C is the radiation pressure cross section 3x3 matrix, unit vector e'l is directed along 
the path of the incident radiation (it is supposed that beam of photons propagate in 
parallel lines) and its orientation corresponds to the orientation of light propagation; F'^ 
is emission component of the radiation force acting on the particle (see Eq. (6)). 

Stationary frame of reference 

Our aim is to derive equation of motion for the particle in the stationary frame of refer- 
ence. 

Covariant equation of motion - first case 

Let the components of the pressure cross section 3x3 matrix C be an orthonormal 
basis e^j, e^2> ^63- We may then write 



dE' 



dT 



0, 
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3 




(3) 



where e'l, e'^ and 63 form an orthonormal basis. Similarly 
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C'e', = Y^ (e;iC'ei) e',, . (4) 



fc=i 



On the basis of Eqs. (Al), (A2), (A3) and (A4), one obtains 
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which corresponds Eq. (37). Thus, the covariant form is represented by Eq. (38). 

More straightforward consideration: Let the components of the matrix C be given in 
the orthonormal basis e^, e'2 and 63 - C^; = C" ej, fc, Z = 1 to 3. Then 

C e; = ^ [el C ei) e^. . (6) 

i=i 

Substitutions A' = C e'i,j = 1,2,3 immediately yield Eq. (37) which has already 
been rewritten to the covariant form represented by Eq. (38). 



Covariant equation of motion - second case 

We want to derive an equation of motion for the particle in the frame of reference in which 
particle moves with actual velocity v. We will use the fact that we know this equation in 
the proper frame of reference - see Eq. (Al). 

Let us have a four- vector = {A^,A), where A^ is its time component and A 
is its spatial component. Since generalized special Lorentz transformations do not form 
a group (in general, composition of two generalized special Lorentz transformations is 
not a generalized special Lorentz transformation) . we will consider more general Lorentz 
transformation. This Lorentz transformation can be written as 

a'" = A^'^A'', (7) 

where summation over repeated indices is supposed (and also in all the following equa- 
tions) - 0, 1, 2, 3 for Greek letters and 1, 2, 3 for Latin letters - and its inverse 

A** = A^" A" . (8) 

Important property of the Lorentz transformation 
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X' = AA (9) 
is, that it can be composed of the following two transformations: 
X = LA, 

X' =RA' , (10) 

where L corresponds to generalized special Lorentz transformation and R represents 
rotation in 3-dimensional space: 

a'" = L"^ A" , (11) 
and (see Eq. (13)) 
L\ = l, 

L\ = L\ = -j{v/c)i , i = 1,2,3, 

U ^ = 6,, + (7 - 1) (vh (v), /v\ i = 1,2,3, j = 1,2,3 , (12) 
where 6ij = 1 if i = j and 6ij = ii i ^ j, 

A^" = L^^" X'^ , (13) 
where (see Eq. (14)) 

'-'a — 'let p II 7 ' 

r}ai3 = diag{+ 1,-1,- I,- I) (14) 

and 

L,' =1, 

Lo' = i/ = 7 Wc)i , i = 1,2,3, 

Li ' = kj + (7 - 1) {v)i {v)j /v^ , i = 1,2,3 , i = 1,2,3 , (15) 

R\ = R\ = Q, z = 1,2,3, 

= , z = 1,2,3, J = 1,2,3, (16) 

where Vij can be expressed in terms of Euler angles, and, moreover, orthogonality condi- 
tions are fulfilled: 

3 

^ ^ ^ij ^ik = 5jk , J = 1, 2, 3 , fc = 1, 2, 3 , 
1=1 

3 

Y^Uj Tkj = 5ik , i = 1,2,3, fc = 1,2,3. (17) 
On the basis of Eqs. (A9) and (AlO) we can write 
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K\ = R>^^L\. (18) 
Eqs. (A12), (A16) and (A18) yield 

K\ = L\ = -j{v/c)i , i = 1,2,3. (19) 
Finally, requirement A^^A^j, = A ^^A'^ yields, 

A^A/ = <5^ (20) 

where 5$ = I \i q = v and 5^ = Q ii q ^ v] Eq. (A7) was also used: A'** = K^",^ A"" , 
— A^'^ A^ . 

Incoming radiation 

Applying Eqs. (A7) and (A19) to quantity {Ei/c,Pi) (four- momentum per imit time - 
proper time is a scalar quantity) and taking into account also the fact that Pi = Ei/c ei, 
we can write 

E'i = E,wi, (21) 
where 

= 7 (1 - V ei/c) . (22) 
Using the fact that = {h u, h v ei) for photons, we have 
u' = uwi . (23) 

We have four- vector p^ = {Ei/c,pi) = (l,ei)£'i/c = {l/wi,ei/wi) w\ Ei/c = 
bi wi Ei/c. We have found a new four- vector b'^, which is given as b^ = A^^ hi in 
the proper frame of reference of the particle: bf = (1, e'^). The transformation of space 
components between bi and bi' corresponds to aberration of light. 

For monochromatic radiation the flux density of radiation energy becomes 

S' = n' hv' c; S = nhuc, (24) 

where n and n' arc concentrations of photons (photon number densities) in the corre- 
sponding reference frames. We also have continuity equation 

df.j'' = 0, f = {cn,cnex), (25) 
with current density j^. Application of Eqs. (A7) and (A19) then yields 
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n' = w\n . (26) 
Using Eqs. (A23), (A24) and (A26) we finally obtain 

S' = wlS . (27) 
Eqs. (A21) and (A27) then together give Ei = wi S A', pi=wi S A' ei/c. 

Covariant equation of motion 

Inspiration comes from the fact that space components of four-momentum are written 
as a product with unit vector e[ . We know that this unit vector can be generalized to a 
four-vector 

6^ = (1 / wi,ei I wi) 

wi = 7 (1 - i; • ei/c) . (28) 

Moreover, we know that w\ S / c is a scalar quantity - invariant of the Lorentz transfor- 
mation (see Eqs. (A27)). 

The idea is to write covariant equation of motion in the form 

^ = 6, . + i t F;. (c - u^) , (29) 

if the result for the emission component of the radiation force acting on the particle was 
added (sec sections 2 and 3 in the main text). The only problem is to find components 
of the tensor of the second rank 

In order to find , we will proceed in two steps. At first, we will rewrite Eq. (A29) 
in the proper frame of reference of the particle. Comparison with Eq. (Al) will yield 
components of G The second step is transformation from G to G^ 

In the proper frame of reference, Eq. (A29) yields 

3 



= ^ }g'^° - Vg'o^ (el), 
^ = ^Ig'^^ -±G'^^ (ei) . \ + ±K, (e;.), , (30) 



i=i ) i=i 

where the term 1 /w'l in brackets is omitted due to the simple fact that it equals 1 in the 
proper frame of reference. Comparison with Eq. (Al) yields 

G' 00 ^q' 3 ^q' j ^ j = 1, 2, 3 , 

G'''^ = -C'kj, j = 1,2,3, fc = 1,2,3. (31) 
In order to find G^ " ,we have to use the Lorentz transformation 
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" = A„'' A^'' g' " , (32) 
where 

■'*-a — '/a p '/ 7 ' 

r]a = diag{+ 1, - 1, - 1, - 1) . (33) 

Usage of the generahzed special Lorentz transformation (Eqs. (A12) or (A15)) yields 
the results of Kimura et al. (2002) - the authors take the space tensor C as a scalar quan- 
tity under Lorentz transformation. However, generalized special Lorentz transformation 
has not to be used. The reason is that any body in torque-frce, accelerated motion un- 
dergoes rotation. This effect is known as the Thomas precession (sec c. g., Robertson and 
Noonan 1968, pp. 66 69). Thus, the process of complete rclativistic derivation and the 
corresponding result presented by Kimura et al. (2002) is incorrect. We can mention that 
the last statement is evident already at first glance: Kimura et al. (2002) take the radi- 
ation pressure cross section 3x3 matrix - space tensor - as a relativistically invariant 
quantity. 

Consistency of the covariant formulations 

We have obtained equation of motion in the form of Eq. (A29). Another form of covariant 
equation is presented in Eq. (35) (Eq. (38)). Are these equations consistent? 

Eq. (38) is covariant equation of motion in the form corresponding to {f3^ = /c) 



7) =^ E^H^^ - + it^eA'^b^-') ■ (34) 

This appendix has discussed an equation of motion of the form 

^G^'^h. + it^e^i^b^-un • (35) 



d 



a T J T, c c , 



We want to show that Eqs. (A34) and (A35) are equivalent, i. e., that {dp^ / d t)j = 

{dp'^/dT)jj. 



We will not write the terms F^, in what follows. 

Multiplication of Eq. (A34) by four-vector bk ^ (and summation over /i) yields 

^) bk^ = ^^^^^j2Q'Abjbk^-l3^bk^),k = l,2,3. (36) 

In calculations of bk f_i we will use the fact that it represents scalar product of two 
four- vectors. Thus, its value is independent on the frame of reference. For the proper 
frame of reference 
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/3'^6fe^ = l , fc = 1,2,3. (37) 
It can be easily verified, that 

6^^ &fe ^ = 1 - e; • e'fe = 1 - (5,- ft , J = 1,2,3 , fc = 1,2,3 , (38) 

since in the optics of scattering processes it is assumed (defined) that unit vectors e'^ , e'2 
and 63 are orthogonal. Thus, we obtain (inserting results of Eqs. (A37) and (A38) into 
Eq. (A36)) 

dp^\ ^ fc = 1,2,3. (39) 



dr J J 

Multiplication of Eq. (A35) by four-vector bk ^ (and summation over /x) yields 

^\ b, =!dlG^-bi,bk^,, fc = 1,2,3. (40) 
dT-Jii c 

Again, the value is independent on the frame of reference. For the proper frame of refer- 
ence Eq. (A31) yields 

,2 



(,77;,/'='^= -^C,.,(ei).(eft), 



(e'fe)^ (C e'l) , fc = l,2,3. (41) 

It was already shown (sec Eq. (A6) and the text below it) that right-hand sides of Eqs. 
(A39) and (A41) are identical. Thus, also left-hand sides of Eqs. (A39) and (A40) are 
identical: 

If we take into account that four-vectors 6^ may be taken in various ways, Eq. (A42) 
yields 

Thus, Eqs. (A34) and (A35) are equivalent, Q. E. D.. 
Appendix B: Einstein's example 

(Reference to equation of number (j) of this appendix is denoted as Eq. (B j). Reference 
to equation of number (i) of the main text is denoted as Eq. (i).) 

Let us consider a plane mirror moving (at a given moment) along x-axis (system 
S) with velocity v = (w,0, 0), v > 0; the mirror is perpendicular to the x-axis (the 
plane of the mirror is parallel to the yz-plane). A beam of incident (hitting) photons is 



53 



characterized by unit vector S' = (cos 9', sin 9', 0) in the proper frame (primed quantities) 
of the mirror. Reflected beam is described by the unit vector e' = (— cos^', sin^', 0) (in 
the proper frame S'). 

The problem is: Find equation of motion of the mirror in the frame of reference S. 

Solution 1: trivial manner 

Consider one photon (frequency /') in the proper frame of the mirror. Since the directions 
(and orientations) of the incident and outgoing photons are characterized by 

S' = (+ cos6»',sin6i',0) , 

e' = (- cos 0', sin 6*', 0) , (1) 
we can immediately write 

= — (!,+ cos sin ^',0) , 



= -L- {\,- COS 0', sin 6*', 0) , (2) 

for the four-momentum of the photon before interaction with the mirror and after the 
interaction. 

As a consequence, the mirror obtains four-momentum 

p'" = Vt - Po=^ 2 cos 9', 0, 0) . (3) 
Application of the special Lorentz transformation to Eq. (B3) yields 
= ^ 27 (cos ^') 1,0,0) , (4) 

where, as standardly abbreviated, 

7 = 1/ ^/^^^^ ; p = v / c. (5) 
On the basis of Eq. (B4), we can immediately write equation of motion of the mirror 

^ = ^27 (cos ^0 (A 1.0,0) , (6) 
or c 

where E[ is the total energy (per unit time) of the incident radiation measured in the 
proper frame of reference. 
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Solution 2: application of general theory presented in Sec. 3 of the main text 

We have to choose orthonormal vectors in the systems S': we will use e'^ = S' and e'2 
and one can easily find 

e'l = S" = (+ cos 61', sin 6*', 0) , 

62 = (- sin 6*', cos 61', 0) . (7) 
We have to write (Qg = 0) 

p' = ^ {Q\ el + Q'^ e^) . (8) 

On the basis of Eqs. (B3), (B7) and (B8) we have (Q(j = 0) 

Q'l = 2 (cos 61')^ , Q2 = - 2 (sin 61') (cos 61') . (9) 
Other prescription yields (see Eq. (34)) 

6? = 7 (l + t;.e'i/c) =7(l + /3cos6i') , 

61 = e'l + [(7 - l)v ■ e'l/v^ +7/c] v = (7 cos 6*' + 7/?, sin 6*', 0) , (10) 
hi = 7 (1 + • e^/c) = 7(1-/3 sin ^'), 

62 = 62+ [(7-l)t;-e2/'u^ + 7/c]t; = (-7sin6'' + 7/3,cos6i',0) . (11) 

Inserting Eqs. (B9) - (Bll) (and = for j = 1, 2, 3) into Eqs. (29) - (30), one obtains 

^ = ^ {[2{cose'f] (65'-/3^) + [-2(sin0')(cos0')] (6^-/3^)} 

= ^27 (cos e')(/3, 1,0,0) . (12) 

Unit vectors e'^ = S' and are used. They are orthonormal in the system S'. However, 
corresponding vectors are not orthogonal in the system S. 

ei = ^ {e'l + [(7 - 1) w • e'l/v'^ + 7/c] v] , 

w 

w' = 7(l + t)-e'i/c) , (13) 

and analogous equation holds for vector 62. Inserting Eqs. (B7), one obtains: 

_ f cos6'' + /3 sin 61' 

~ tl + /3cos6l''7(l + /3cos6")' 
f -sing' + /3 COS0' 



l-^sine*' '7(l-/3sin0') 
It can be easily verified that scalar product of these two vectors is nonzero, in general. 
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Comparison with Einstein's result 

Inserting = S A'^i^ror cos6»' into Eq. (B6) (or Eqs. (B12)) and using Eq. (B14) 
(d = (cos^, sin^, 0)) for the purpose of obtaining cos^' = (cos^ — /3)/(l — /3 cos^), one 
easily obtains: 

i) dE/dr = 2 7^ 5 ^^irror (cos^ — (3)^ (3; using definition of radiation pressure dE/dt = 
P V A'^irror, we have P = 2 (5 / c) (cos^ - / (1 - p^), or, 

ii) dp/dr = 27^ (5 ^mj^roT- / ^) (cos^ — P)^', using definition of radiation pressure P = 
[dp/dt) I A'^,„^^, we have P = 2{S /c) {cosO - / (1 - /J^). 

Result for P is consistent with the result presented in Einstein (1905). 

Appendix C: Equation of motion presented by Kimura et al. (2002) 

(Reference to equation of number (i) of this appendix is denoted as Eq. (C i). Reference 
to equation of number (j) of the appendix B is denoted as Eq. (B j). Reference to equation 
of number (k) of the main text is denoted as Eq. (k).) 

It is said that the paper by Kimura et al. (2002) is relevant; irrelevant are papers by 
Klacka (2000a, 2000b, 2000c), or derivations and results presented in Sec. 3 of this paper. 
Thus, it is important to clariiy the situation by presentation of detailed arguments. 

Kimura et al. (2002) derive and present equation of motion of the form 

dv S A' 
dt m c 
C = Q[ [{I - V- ki/c) fei - v/c] + 

Q'^lil - 2v ki/c) k2 + {v- k2/c) fei] + 

Q'3 (1 - 2vki/c) fca , 

e[ = (1 + V ■ k\/c) k\ — v/c, 

e'2 = k2 + {v ■ k2/c) ki , 

e'3 = ka , 

e» • e'j = S,j , i,j e {1,2,3} , 

ki • kj = Sij , i,j G {1,2,3} , 

vks = 0, (1) 

(see Eqs. (1), (2), (3), (10), (14) in Kimura et al. 2002) and it is supposed that F'^ = 0. 
There have appeared two independent suggestions: 

i) Eq. (CI) is correct and Eq. (36) is incorrect. 

ii) Eq. (CI) is equivalent to Eq. (36). As a proof, the following argument is presented. 



56 



Equations of motion presented by Klacka (2000a, 2000b, 2000c) and Eq. (36) (for the 
case = 0) and Kimura et al. (2002) arc equivalent and all the difference conies from 

usage of different basic vectors. The set of non-orthogonal unit vectors {ej^j = 1,2,3} 
(see Eq. (36)) is replaced by new set of orthonormal vectors {kjij = 1,2,3} in the way 

ei = fei , 

62 = k2 + iv ■ k2/c) (fei - ^2) + V/C , 

63 = k3- {v ■ ka/c) k3 + v/c , 



Eq. (C2) yields e[ = {1 + v ■ ki/c)k\ — v/c, e'2 = k2 + {v ■ k2/c)ki and 63 = ks, which 
is equivalent to Eq. (CI). 

Kimura et al. (2002) and Einstein's example 

Let us consider a plane mirror moving (at a given moment) along x-axis (system S) with 
velocity v = (w,0, 0), w > 0; the mirror is perpendicular to the x-axis (the plane of the 
mirror is parallel to the yz-plane). A beam of incident (hitting) photons is characterized by 
unit vector S' = (cos sin ^',0) in the proper frame (primed quantities) of the mirror. 
Reflected beam is described by the unit vector e' = (— cos^', sin^', 0) (in the proper 
frame S'). 

We are interested in application of Kimura et al. (2002) general equation of motion 
to this example. 

On the basis of considerations presented in Appendix B, Eqs. (B7) and (B9) imme- 
diately yield 



V k^ = . 



(2) 



= (+ cos6'',sin6»',0) , 



62 = (- sin 6*', cos 6*', 0) . 

Q'l = 2 (cos 61')^ , <52 = - 2 (sin 6'') (cos6i') , Q3 = . 



(3) 



Eq. (CI) yields for orthonormal vectors 



fei = (1 — v-e'i/c)e'i + v/c, 



k2 = e'2 - {v ■ e'a/c) e[ . 



(4) 



Using Eqs. (C3) and (C4) 




(5) 
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Inserting Eqs. (C3) and (C5) into Eqs. (CI), one easily obtains 

1 + 2 (cos 6*')^ ,0,o| . (6) 



djv_ ^ SA'^,,,^, cos 9' r _ v 

at TO c L c 



Using the relation = S" (1 - 2 v cos9'/c) A'^^^^^^ cos 6*' in Eqs. (B6) or (B12), 
the correct result can be written in the form 



|2 (l - 2 ^ cos 6*') cos6i',0,o| . (7) 



rft; ^ S A!^„^^^ cose' 
dt m c 

Eq. (C6) was obtained on the basis of equation of motion presented by Kimura et al. 
(2002), while the correct result corresponds to Eq. (C7). Thus, Kimura et al. (2002) are 
incorrect - their general equation of motion yields result not consistent with Einstein's 
result (Einstein 1905). 

Physics of transformation represented by Eq. (C2) 

Transformation {e[\i = 1,2,3} — > {ei;i = 1,2,3} from frame of reference S' to frame 

of reference S corresponds to aberration of light (see Eq. (36) or heuristic derivation 
Sec. 3.4). In the given frame of reference S we may want to define two sets of unit 
vectors {e^; i = 1, 2, 3} and {U; i = 1, 2, 3}, where ei = li and • ej ^ Sij, It ■ Ij — Sij, 
i,j G {1,2,3}. Relation between = 1,2,3} and {ei;i = 1,2,3} corresponds to pure 
spatial rotation (geometry) and it can be described as 

ei = h , 

62 = £,ili + 6 h + 6 ^3 , 

63 = ??1 + 112 h + ??3 ^3 • (8) 

The conditions 

h-lj = Sij , i,ie {1,2,3} (9) 
lead to 

^1 = ei • 62 , ryi = ei • 63 , 

3 

E = 62 • 63 . (10) 

The important property is that both sets of unit vectors {li;i = 1,2,3} and 
{ef, i = 1,2, 3} are defined in the same frame frame of reference S. As for transformation 
defined by Eq. (C2), the sets of unit vectors {6,; i = 1,2, 3} and {fej; i = 1,2, 3} are also 
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defined in the system S. But what about translational velocity terms + v/cin Eq. (C2)? 
Any velocity term + v/c in transformation corresponds to transformation from one (lo- 
cal) incrtial frame of reference to another one. Thus, the conclusion is evident: Eq. (C2) 
is of no physical sense. 
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